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Abstract 

In this paper, we first study a new sensitivity index that is based on higher moments and gener¬ 
alizes the so-called Sobol one. Further, following an idea of Borgonovo (ED, we define and study a 
new sensitivity index based on the Cramer von Mises distance. This new index appears to be more 
general than the Sobol one as it takes into account, not only the variance, but the whole distribution 
of the random variable. Furthermore, we study the statistical properties of a Monte Carlo estimate 
of this new index. 

Keywords: Sensitivity analysis, Cramer von Mises distance, Pick and Freeze method, func¬ 
tional delta-method, Anderson-Darling statistic. 


1 Introduction 


A very classical problem in the study of computer code experiments (see |2T]) is the evaluation of the 
relative influence of the input variables on some numerical result obtained by a computer code. This 
study is usually called sensitivity analysis in this paradigm and has been widely assessed (see for example 
m, m, E] and references therein). More precisely, the result of the numerical code Y is seen as 
a function of the vector of the distributed input {Xr)r=i,--- ,d {d S N*). Statistically speaking, we are 
dealing here with the unnoisy non parametric model 

y = /(Ai,...,Vd), (1) 

where / is a regular unknown numerical function on the state space Ei x E 2 x ... x E^ on which the 
distributed variables (Vi,..., Xfi) are living. Generally, the random inputs are assumed to be independent 
and a sensitivity analysis is performed by using the so-called Hoeffding decomposition (see [23] and [T]). 
In this functional decomposition, / is expanded as an L^-sum of uncorrelated functions involving only a 
part of the random inputs. For any subset v of Id = {1,..., d}, this leads to an index called the Sobol 
index (\22\) that measures the amount of randomness of Y carried in the subset of input variables (Xi)i^y. 
Since nothing has been assumed on the nature of the inputs, one can consider the vector {Xi)i^y as a 
single input. Thus without loss of generality, let us consider the case where v reduces to a singleton. The 
numerator Hy of the Sobol index related to the input Xy is 


Hy 


Var (E [V|V„]) = Var(V) - E 




E[F|V„])" 


( 2 ) 


while the denominator of the index is nothing more than the variance of Y. In order to estimate Hy 
the clever trick discovered by Sobol [32] is to rewrite the variance of the conditional expectation as a 
covariance. Further, a well tailored design of experiment called the Pick and Freeze scheme is considered 
m- More precisely, let X"" be the random vector such that A)) = Xy and A" = A' if i u where A' is 
an independent copy of A^. Then, setting 


A" := /(AG (3) 

an obvious computation leads to the nice relationship 

Var(E(y |A„)) = Cov (A, A"). (4) 
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The last equality leads to a natural Monte Carlo estimator (Pick and Freeze estimator) 


rjiV _ 


1 ^ 

i=i 


N 


2N 




i=i 



( 5 ) 


where for j = 1, ■■■, N, Yj (resp. Y^) are independent copies of Y (resp. F"). The sharp statistical 
properties and some functional extensions of the Pick and Freeze method are considered in Hg, dS] and 
m- Notice that the Sobol indices and their Monte Carlo estimation are based on order two methods 
since they derived from the Hoeffding functional decomposition. This is the main drawback of this 
kind of methods. As an illustration consider the following example. Let Xi and X 2 be two independent 
random variables having the same first four moments (equal e.g. to 1) and such that E [Af] ^ E [Xf]. 
Let us consider the following model 

r = Ai + A2 + xfxl 


Then 

Var (E [F|Ai]) = Var(Ai + Xf) = Var(A 2 + A|) = Var (E [Y\X 2 ]) . 

However, since A is a symmetric function of the inputs Xi and X 2 that do not share the same distribu¬ 
tion, Xi and X 2 should not have the same importance. That shows the need to introduce a sensitivity 
index that takes into account not only the second order behaviors but all the distributions. As pointed 
out before, Sobol indices are based on decomposition. As a matter of fact, Sobol indices are well 
adapted to measure the contribution of an input on the deviation around the mean of Y. However, it 
seems very intuitive that the sensitivity of an extreme quantile of Y could depend on sets of variables 
that cannot be read only in the variances. Thus the same index should not be used for any task and we 
need to define more adapted indices. There are several ways to generalize the Sobol indices. One can, 
for example, define new indices through contrast functions based on the quantity of interest (see HSI). 
Unfortunately the Monte Carlo estimator of these new indices are computationally very expensive. In 
[ 9 ], the author presents a way to define moment independent measures through dissimilarity distances. 
These measures define a unified framework that encompasses some already known sensitivity indices. 
Unfortunately, the estimation of such indices relies on the estimation of density ratio estimation that 
can be computationally expensive. Now, as pointed out in 0,0, m, CH] and m, there are situations 
where higher order methods give a sharper analysis on the relative influence of the input and allow finer 
screening procedures. Borgonovo et al. propose and study an index based on the total variation distance 
(see 0, m and Cl)- While Owen et al. suggest to use procedures based on higher moments (see |18| . 
m)- Our paper follows these tracks. We will first revisit the works of Owen et al. by studying the 
asymptotic properties of the multiple Pick and Freeze scheme proposed therein for the estimation of 
higher order Sobol indices. Further, we propose a new natural index based on the Cramer von Mises 
distance between the distribution of the output Y and its conditional law when an input is fixed. We 
will show that this approach leads to natural self-normalized indices as in the case of the Sobol-Hoeffding 
decomposition of the variance. As a matter of fact, as for Sobol indices, the sum of all first order indices 
can not exceeds one. Notice that these indices extend naturally to multivariate outputs. Furthermore, 
we show that surprisingly a Pick and Freeze scheme is also available to estimate this new index. The 
sample size required to build such an estimator is of the same order as the size needed for the classical 
Sobol index estimation allowing its use in concrete situations. 


The paper is divided in three sections. In the next section, we will study the statistical properties of the 
multiple Pick and Freeze method proposed earlier by Owen et al (HE], HSj). Section is devoted to the 
new index built on the Cramer von Mises distance. In the last section, we give some numerical simulation 
that illustrate the interest of the new index. In particular, we revisit a real data example introduced in 
0 and studied in [TE] and 0. 


2 Multiple Pick and Freeze method 


Using the classical Hoeffding decomposition, for a singleton v G Id, the numerator of the classical Sobol 
index with respect to v is given by 


Hf=E 


(E[r|A„] - E[y])' 


( 6 ) 
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Following [T^ and m, we generalize this quantity by considering higher order moments. Indeed, for any 
integer p ^ 2, we set 


(7) 


( 8 ) 


ff^=E[(E[r|x„]-E[r]f]. 

Hy = Hy. The following lemma gives the Pick and Freeze representation of for p ^ 2. 

Lemma 2.1. For any v € Id, one has 

p 

E [(E[F|X„] - E[y])P] = E (Y^’^ - E[r]) 

j=i 

Here, = Y and for i = 2,... ,p, P’'’* is constructed independently as Y^ defined in equation (|^. 
Obviously, Hf is non negative for even p and 

\HP\^E[\Y -E[Y]f>]. 

Further, is invariant by any translation of the output. 


Estimation procedure In view of the estimation of Hf, we first expand the product in the right-hand 
side of (Isl to get that 


1^0 ^ ^ 


nr- 


with the usual convention = 1- Second, we use a Monte Carlo scheme and consider the following 

Pick and Freeze design constituted by the following p x A^-sample 

(F-). 


ii,j)eipXiN 


We define for any any N G N*, j G In and I G Ip, 


pv _ P 


-1 


I 


1 


N 


ki<...<ki£lp \i—l 

The Monte Carlo estimator is then 




TTV _ 

^p,N — 


t(^)(-ir-(pr) 


_oAP-^—7; 


Pi- 


(9) 


Notice that we generalize the estimation procedure of m and use all the available information by 
considering the means over the set of indices ki,..., ki G Id, kn k^n- The following theorem provides 
asymptotic properties of 

Theorem 2.2. ^ is consistent and asymptotically Gaussian: 

Vn{hIn-h;) 


( 10 ) 


where 


a^=p [Var(y) + {p - l)Cov{Y, W’^)] ( ^ aik ) , 
I 


,1=1 


ai = -E[Yy-fi l = l,...,p 

P 

p-i 


and 


6i = (-l)^-V(p - l)E[Ff-i -I- ^ {-ir-‘ip - /)E[r]^'-'-iE 

1^2 ^ ^ 

Q(-i)p-'E[r]p-', i = \,...,p. 




.1 = 1 


bi = 
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Proof of Theorem \2.S\ The consistency follows from a straightforward application of the strong law of 
large numbers. The asymptotic normality is derived by two successive applications of the delta method 

m ■ 


(1) Let IT/ = (F/ 
is given by 


Yj^pyr (j = 1,..., N) and the mapping from to whose l-th coordinate 


.. ,Xp) 



E 

ki < ... <ki 


ki € Ip , Z 1J . . . , / 



Let be the covariance matrix of Wj. Clearly, one has Ef = Var(F) for i G Ip while EL = 
Cov(y*'’®, = Cov(y, The multidimensional central limit theorem gives with m = (E[y],..., E[F])^ 


N 


We then apply the so-called delta method to and g^ so that 

(w]v) -g^ (E [W^])) ^4^Af(0,Jgi (E [W^])Y^Jy (E [W^]f 


with Jgi (E [VL^]) the Jacobian of g^ at point E [W^]. Notice that for i G Ip and k G Ip, 

(E [W^]) = = -E[Yf-^ =: ai. 

OXk ^ ^ (/) p 

Thus E^ := Jgi (E [W^]) EV^i (E [W^])^ is given by 


Yjg = pa,aj (E}i + {p - . 

(2) Now consider W/ = (j = 1,..., N) and g^ the mapping from to K defined by 

We apply once again the delta method to so that 

^ {y^ {wl) - 9^ (E [W2] )) ^4^ M (o, Jg. (E [W2] ) EVg. (E [W2] )^) 


with Jg 2 (^E [JF^]) the Jacobian of g^ at point E [JT^]. Notice that for k G Ip, 

(E [w^]) = (-if-V4-i)E[rf-i 


dg^ 


dyi 


+ Y(f^i-lY-‘(p-l)E\Yr-'-'E 


n>'- 


and 


dyi 


(E [W^]) = Q(-lf-'E[r]P-'. 

Thus the limiting variance is 

4 := Jg2 (E [W^]) Y^Jg. (E [W^])^ =p{Y\^ + (p- 1)E}2) , 

where bi is the i-th coordinate of V g^ (E [W'^] ). 


\i=l 


□ 
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The collection of all indices Rp is much more informative than the classical Sobol index. Nevertheless 
it has several drawbacks: it may be negative when p is odd. To overcome this fact, we may have 
introduced E[|E[F|Xi,z G ?;] — E[y]|^] but proceeding in such a way, we would have loose the Pick and 
Freeze estimation procedure. The Pick and Freeze estimation procedure is computationally expensive: it 
requires a p x N sample of the output Y. In a sense, if we want to have a good idea of the influence of 
an input on the law of the output, we need to estimate the first d indices and hence we need to run 
the black-box code K x N times. Moreover, these indices are moment based and it is well known that 
they are not stable when the moment order increases. In the next section, we introduce a new sensitivity 
index that is based on the conditional distribution of the output and requires only 3 x N. 


3 The Cramer von Mises index 

In this section the code will be denoted by Z = f{Xi ,..., Xii) G Let F be the distribution function 
of Z. For any t = (ti,... ,tk) G 


and F'"(t) the conditional distribution function of Z conditionally on 

F^(t) = P (Z ^ ) = E [1 |V„] . 

Notice that {Z ^ t} means that {Zi ^ t\^...,Zk ^ tk}- Obviously, E[F"(t)] = F(t). Now, we apply 
the framework presented in Section^ with Y{t) = M and p = 2. Hence, for t G fixed, we have a 

consistent and asymptotically normal estimation procedure for the estimation of 


E 


{F{t) 


F^{t)f 


We define a Cramer Von Mises type distance of order 2 between C (Z) and £ (ZjXy) by 


d: 


2,CVM 


:= / E 


{F{t)-F^{t)f dF{t). 


( 11 ) 


The aim of the rest of the section is dedicated to the estimation of cvm study of the asymptotic 

properties of the estimator. Notice that 


D 


V 

2,CVM 


E 


E 


{F{Z)-F^iZ)y 


( 12 ) 


Let us note that these indices are naturally adapted to multivariate outputs. 

Remark 3.1. Lfnlike the procedure for p = 2, we did not normalize the generalized Sobol index of Y{t). 
The purpose, that becomes clear in this section, is to avoid numerical explosion during the estimation 
procedure. Indeed, the normalizing term would be F{t){l — F{t)), like in the ANderson-Darling statistic, 
canceling for small and large values of t. Nevertheless, in view of the following proposition, one can 
consider qvm instead of D 2 cvm i'^ order to have an index bounded by 1 as for the Sobol index. 
The asymptotic properties will not be affected by this renormalizing factor, so we still consider D 2 cvm- 

Proposition 3.2. One has the following properties. 

1. 0 ^ F) 2 CVM ^ Moreover, if k = 1 and F is continuous, we have 0 ^ FI^cvm ^ 

2. D 2 CVM invariant by translation, by left-composition by any nonzero scaling ofY. 

We then proceed to a double Monte-Carlo scheme for the estimation of cvm consider the following 
design of experiment consisting in: 

1. two fV-samples of Z: 1 ^ j ^ TV; 

2. a third TV-sample of Z independent of {Zj'^ , Wk, 1 ^ fc ^ TV. 
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The empirical estimator of qvm then given by 


D 


V 

2,CVM 


N 


N 

E 


1 

N 


N 

{Z]-^siWk} 

i=i 


1 

m 


E( 

t=i 


{Z]’^!iWk} + ^{Z]'^!iWk} 


The consistency of qvm follows directly from the following lemma: 


Lemma 3.3. Let G and H be two — measurable functions. Let (Uj)j^i^ and (T4)fcg/„ be two inde¬ 
pendent samples of iid rv such that E[G(C/i, Vi)] = 0 and E,[H(Ui,U 2 ,Vi)] = 0. We define Sn and Tn 
by 


1 - 1 ^ 

E ^ E H{U^,Uj,Vk). 

j,k — l i,j,k—l 


Then Sm and Tjv converge a.s. to 0 as N goes to infinity. 


Proof, (i) If we prove that E[5'^] = O (^), we then apply Borel-Cantelli lemma to deduce the almost 
sure convergence of to 0. Clearly, 




where the sum is taken over all the indices ii, i 2 , is, i^, ji, j 2 , js, Ja from 1 to N. The only scenarii that 
could lead to terms in O (;^) or even O (1) appear when we sum over indices all different except 2 i’s or 
2 j’s or over indices all different. Nevertheless, in those cases, at least one term of the form E[G(C/i, Vj)] 
appears. Since the function G is centered, those scenarii are then discarded. 

(ii) Analogously, it suffices to show that E[T^] = O (^). The only scenarii that could lead to terms 
in O (^) or even O (1) appear when we sum over indices all different except 2 i’s, 2 j’s or 2 k’s or over 
indices all different. Nevertheless, in those cases, at least one term of the form E[H{Ui, Uj, 14)] appears. 
Since the function H is centered, those scenarii are then discarded. □ 


Corollary 3.4. cvm strongly consistent as N goes to infinity. 


Proof. The proof is based on Lemma 

F{Zj,Wk) = I fo ^{zj 


First, we define , G{Zj,Wk) = 

and H{Z,,Z^,Wk) = F{Z,,Wk)F{Z^,Wu). 


^ {Zf^siWk}^ {Z]'^!iWk}’ 
Second we pro- 
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ceed to the following decomposition 


d: 


2,CVM 


1 ^ 




1 ^ 

3 = 1 


1 ^ 

i=i 




>) 


N 


= E ", 


N 

E ( 


7V2 ^ {z]-^^Wk} 4jY3 

AT AT 

= ^ E E HiZ,,Z„W,) 


II 




)( 


II, 


+ 1. 


j,k=l 


i,j,k—l 


1 1 

= ]^ E H^fe) - EIG(Z„ Wk)]} - E M"fc)]} 

TV TV 

+ ]^ E E E[i/(Z„Z,,ty,)] 


i.fc=i 


2 ,ji,fc = l 


1 i\ 1 

= E {GiZ^Wk) - E[G(Z„ Wk)]} - E W^fc)]} 

j,k—l i,j,k—l 

+ E[G{Z^,Wi)] - (l - ^) miZi, Z2,Wi)] - ^E[H{Z,, Z,,Wi)]. 


The two first sums converges almost surely to 0 by Lemma 3.3 The remaining term goes to E[G{Zi^ hhi)] — 
E[H{Zi^ ^ 2 , Wi)] as N goes to infinity. 


It remains to show that qvm ~ ^[^(^ 1 ) ^i)] ~ ^ 2 , W^i)]- On the one hand, 


DIcvm = [ E[(F(t) - F^{t)f]dFit) = E[H^^{W)] 

JR 

= E[Cov(ll { 2 ”-"s;wi}’ ^ {zf-^sswi})] 

= Eiy[Ez[ll {z”,"^Wi}ll ~ Ez[^ 


On the other hand, 


E[G{Zl,W^)]-E[H{Z^,Z2,W^)] 

= - -E[^ll + M (ll {z”-"^VKi} + ^ {z^’^^Wi})] 

= Ew[E2 [1 I {z”.i^vvj}1 I {2 ”-^<;wi}]] ~ E[ll {z^-^^Wi}] 

= Ew[E2[1I {z”.i^Wj}1I {2”-^^Wi}]] ~ E[E[1I {z”.i^Wj} 1I {2j'^^w'i}l^i]] 

= Ew[E2[1I {zi'^t^Wi}^ {zl''^t^Wi}\\ ~ E[E[1I ^2”,"^vVi}I^i]E[ 1I {zj’^^iVi}l^i]] 

= Evi/[E2[1I {zi’^^Wi}^ {zi'^^Wi}^ — E[E[1I ^2”,"^vvj}|W^i]]E[E[ll 
= Evi/[Ez[ll {z^.i^vvi}^ ~ {zL^^rvi}]E[ll {zj-^^Wi}] 

= Evi/[Ez[ll {z”.i^ivj}ll {2”.2 ^vVi}]] ~ E[ll 
= Ew[E2[1I {z”.i^Wj}1I {2”-^^Wi}] ~ E2[1I {2”’"^Wi}]^]' 


□ 

We now turn to the asymptotic normality of cvm ■ follow van der Vaart [23] to establish the 

following proposition (more precisely Theorems 20.8 and 20.9, Lemma 20.10 and Example 20.11). 
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Theorem 3.5. The sequence of estimators qvm asymptotically Gaussian in estimating Z ?2 cvm 
that is cvm ~ ^2 cvm'^ weakly convergent to a Gaussian centered variable with variance 

given by ( |l^ . 

Proof. We define 


N 




i=i 

. N 


i = l,2, 


1 ^ 

^N{t) = 


fc=i 


and rewrite qym ^ ^ regular function depending on the four empirical processes defined behind: 


D 


2,CVM 




G]y + 


dF 


N- 


Since these processes are cad-lag functions of bounded variation, we introduce the maps ipi, 4>2 '■ 
i?Vi[— 00 ,-l-oo]^ !-)■ M and : i3Vi [— 00 ,-|-oo]^ >-)• K by 

MFi,F 2 )= JiF^ydF^ and ^(Fi, F 2 , ^^ 4 ) = ^^ 4 ) - ^^2 ^^ 4 ) , 

where set BVmIo, b] is the set of cad-lag functions of variation bounded by M. 

By Donsker’s theorem, 


v^(Gy-F,G%-F,G]f-G,¥N-F) 4 G 

where G{t, s) = P ^ t, Z'"'^ ^ s), G{t) = G{t, t) and G is a centered Gaussian process of dimension 
4 with covariance function defined for (t, s) G by 

n(<, s) = E (W^J) - E (X,) E 

and Xt := (ll {z»,i^t}, II H {Z’'4^t}ll H {rv^t}) • 

Using the chain rule 20.9 and Lemma 20.10 in [5^, the map is Hadamard-differentiable from the 
domain 00 , -|-oo]^ into K. The derivative is given by 


(hi,h2,h3,h4) H> 4(^3,_F4)(^ 344 ) - 



where the derivative of if (resp. (/>) are given by Lemma 20.10: 




h 2 -d(p o Fi + 


ip'{Fi)hidF2 


taking ip = Id (resp. p{x) = x^) and /i_ is the left-continuous version of a cad-lag function h. 
Since 

F2,cvm — (g]v,G^,G^ jEa?) 5 








we apply the functional delta method 20.8 in [23] to get limit distribution of \/N [D^ cvm ~ ^2 cvm 
converges weakly to the following limit distribution 


h4-d{F^ -G)+ hsdF 


F{hi + h2)dF. 


Since the map is defined and continuous on the whole space 00 , + 00 ]^, the delta method in its 

stronger form 20.8 in implies that the limit variable is the limit in distribution of the sequence 




{F,F,G,F) 

= Vn 


N{Gl,-F,G%- F,G\^^ - G,¥n - F 


)) 


J (Fjv - F)_ d (f^ - G)) + J (G]f -G-F{Gj,+Gjf- 2F)) dF 


We define 


U-.= J )\{w<t}d{F\t)-Git,t))=G{W+,W+)-F{W+f, 

V:= J + ^{z^.2^t}) F^t)] dF{t) = ^ {F{Z-’^f + FiZ^’^) - F{Z'’^^ V 

Obviously, 


E(G) = J (G(4,4)-F(4)2) dF{t), 

E{U^)= I {GiU,U)-F{u f)^dF{t) 

E(F)= J {F{tf - G{tG)) dF{t), 

E(^") = 11 F{tfdF{t) + 11 
By independence, the limiting variance is 

= VarG + Vary. 


F{t V s) {F{t Vs) - F{tf - F{sf) + -F{tfFisf 


dG{t,s). 


(13) 

□ 


4 Numerical applications 

4.1 A flavour of the method in a toy model 

Let us consider the quite simple linear model 

V = crVi + W 2 , o; > 0, 

where Xi has the Bernoulli distribution with success probability p and Xi, X 2 are independent. Assume 
further that A 2 has a continuous distribution F on M with finite variance cr^ and that p = E[A 2 ] and 
cr^ = a^p{l — p). With these choices the random variables aAi and A 2 share the same variances and 
Xi and X 2 have the same first order Sobol indices (1/2). On one hand, the conditional distribution Y 
knowing Xi = 0 is the same as the one of X 2 and the conditional distribution Y knowing Ai = 1 is 
F{- — a). On the other hand, the conditional distribution of Y knowing X 2 is 

P (F = a + A 2 I A 2 ) = 1 - P (V = V 2 IA 2 ) = p. 

Hence, the density of Y is the mixture pF{- — a) + (1 — p)F{-). Tedious computations lead to 

DI,cvm = P(1 -P) [ {F{t) - F{t - a)f [(1 - p)dF{t) + pdF{t - a)] (14) 

jR 
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and 


D2,CVM ~ ~ ~ P) 



a)dF{t) 


(15) 


As p goes to 0 (and a goes to infinity), cvm S^es to 0 and qym goes to 1/6 while the two classical 
Sobol indices remains eqnal to 1/2. Our new indices shed lights on the fact that, for small p, X 2 is much 
more influent on Y than Xi which follows the intuition but is lost when one computes the classical Sobol 
indices. 


Similarly we can compute the indices of order g (g ^ 2): 

= a'^ [p(l - p)"^ + (-p)'^(l - p)] 

Hl=E[{X2-t^n 

Some examples (i) if X 2 is a centered Gaussian with variance cr^ = a^p{l — p), one can easily derive 
an explicit formula for the second index of order g: 


Hi = E[{X2 - my] 


0 

29/2.(q/2)! 


if g is an odd number 
else. 


(ii) if X 2 is a uniformly distributed on [0,6] with b = 2a-\y3p{l — p), one can easily derive an explicit 
formula for the different indices introduced before: 


^2,CVM — P (1 — P) ^ 


(f)'(i-if) 

1/3 else, 


D 


2 

2,CVM 


1 

6 


pj^-p) 

2 




and 


HI = e[(A2 - py] 


0 if g is an odd number 

(6/2)V(<7 + 1) else. 


(iii) if X 2 is a exponentially distributed with mean 1/A = ay/p(l — p), one can easily derive an explicit 
formula for the different indices introduced before: 


ni 

^2,CVM 




y and D. 


2,CVM 




and 


HI = E[(A2 - py] = 

The results are presented in Figures to The blue line (resp. the red dashed line) represents the 
true value of index (resp. P^cvm)- blue line with o (resp. the red dashed line with +) 

represents the estimation of index Qyjy (resp. cvm)- 


4.2 A non linear model 

Let us consider the quite simple model 

Y = exp{Ai + 2 A 2 }, 

where Xi and X 2 are independent standard Gaussian random variables. Straightforwardly, we can derive 
the density function of the output Y and its distribution function: 

/v(y) = -^e-('"*^)^/i°^R+( 2 /) and Fy(y) = $ 

where 4 stands for the distribution function of the standard Gaussian random variable. Its density 
function will be denoted / in the sequel. Then tedious computations lead to the Sobol indices 
and Dl Qyjy. 


10 













0,1 


I 0 05 

C 


0 


J_1_1— 1 —i- 

' ' 1 



—t?--H) 


0 0.5 1 1,5 2 

Value of a for fixed p=1/4 




0 0.5 1 1.5 2 

Value of a for fixed p=1/2 



0,1 


I 0.05 

C 


0 


-I- I I I 

G-^-6-e- 


0 0.5 1 1.5 2 

Value of a for fixed p=3/4 


Figure 1: Example 1 - X 2 Gaussian distributed. 
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Figure 2: Example 1 - X 2 uniformly distributed. 


11 


































0 0.5 1 1.5 2 

Value of a for fixed p=1/4 




0.1 

0.05 





1 1 1 ' * 

r\ 


V -0 




0 0.5 1 1.5 2 

Value of a for fixed p=3/4 


Figure 3: Example 1 - X 2 exponentially distributed. 


Proposition 4.1. 


D. 


1 


1 


2,CVM 


= — arctan2-« 0.019 


and 


D 2 cvM = “ arctan VT9 — « 0.095. 

’ 3 


Proof. First of all, the distribution function of Y\Xi is given by 


FW(t) = P(y ^ t\Xi) = $ 


Then 


D. 


2 ,CVM 


= / E 


{F^^){t)-FY{t)f] fY{t)dt 


= / E 
Vr+ 


= / E 




y/lOirt 






dz 


= E 


$(X2) - $ 


/ 75X2 - 

V 2 - ^ 


where Xi and X 2 are independent standard Gaussian random variables. In the same way, 

DlcvM = '^\iHX2)-^(V5X2-2Xi)f . 


(16) 


(17) 
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Thus we are lead to compute the bivariate function: 


^(a, p) := E [($(^2) - $ {aX2 - PX^f] 
at (a, /3) = (•\/5/2,1/2) and (a, /3) = (-x/S, 2). The term E [$(^ 2 )^] is 


E [$(X2)2] = J ^zpf{z)dz = 




+ 00 


We introduce U, U' and V independent random variables distributed as a standard Gaussian for the two 
first and a centered Gaussian with variance for the third one. Then the term E $ {aX 2 — PXif' 

can be rewritten as 


E 


^{aX2-l3XiY 


= E 


^{vy 


= E 


= E [E [1 u^v\V] E [1 c/-^y|F]] = E [E [1 
= E [1 u'^v] = F{U ^V,U' ^V). 


Let G be the real-valued function defined on R by G{a) = P(C/ ^ aV, U' ^ at^) where U, U' and V are 
independent standard Gaussian random variables. We want to compute G{-\/a‘^ + /3^). Integrating by 
parts, we have 


Jm ) 


Since G(l) = 1/3, we get 


G\a) =2 [ 

Jr 


1 


7r(a^ -I- 1) 

a 


<i)(a2:) 


dz 

2tt 

-(aWl)zV2 


■K{aJ + 1 ) + 1 


G{a) = 


1 

3 


r(a;2 


1 ) 


— dx= —I— (arctan \/l + 2a^ — arctan -s/S) 
1 3 TT 


— arctan \/T+~2aJ 

TT 


and 


E 


$ (aX 2 — fiXxf' = i -I- — (arctan \/l + 2{a'^ + p^) — arctan pS) = — arctan pi + 2{a^ + p^). 
J 3 TT TT 


In the same way, the last term E [<I>(J'f 2 )$ ( 0:^2 — PXi)] is given by 

E [$(X2)$ [aX2 - yXi)] =Fiu^V, \l ^-^U' ^ ^ . 




where U, U' and V are independent standard Gaussian random variables. Remind we only need to 
consider {a^fi) = (•\/5/2,l/2) and {a,P) = {p5,2) in which cases \J= 1. Thus the last equals 1/3 
in both cases that leads to the result. □ 


Remark 4.2. In the previous proof, we show that 


G{a) := P ([/ < aV, U' < aV) = — arctan pl + 2a? 

TT 

where U, U' and V are independent standard Gaussian random variables. Actually, this result is also 
a straightforward consequence of Lemma 4.3 in m at 0 with X = (aV — U) fpaPPT and Y = (aV — 
Uyjpa? Y 1. Nevertheless, since our proof is different and elegant, we decide not to skip it. 
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We can compute the indices of order q {q ^ 2): 


Hf=E 
Hl = E 


(gXi+2 _ g5/2^g 
^g2Xi+l/2 _ g5/2^g' 


The results are in the following tabular. 

Cramer von Mises Sobol indices 


~jTL ff2 ^ 

^2.CVM ^2,CVM ^^ 



^2.CVM 

^2,CVM 

O 

O 

True values 

0.0191 

0.0949 

0.0118 

0.3738 

N = 10^ 

0.0372 

0.0960 

0.1962 

0.1553 

N = 10^ 

0.0192 

0.0929 

0.0952 

0.1085 


As a conclusion, with only N = 10^, the algorithm provides a precise estimation of the different indices. 
Moreover, in this example, Sobol and Cramer von Mises indices give the same influence ranking of the 
two random inputs. Nevertheless, it seems that the estimation of the Cramer von Mises indices is more 
efficient to give the true ranking. 


4.3 Application: The Giant Cell Arthritis Problem 

Context and goal 

In this subsection, we consider the realistic problem of management of suspected giant cell arthritis 
posed by Bunchbinder and Detsky in [S]. More recently, this problem was also studied by Felli and 
Hazen m and Borgonovo et al. [S]. As explained in [H], “ giant cell arthritis (GCA) is a vasculitis of 
unknown etiology that affects large and medium sized vessels and occurs almost exclusively in patients 
50 years or older”. This disease may lead to severe side effects (loss of visual accuity, fever, headache,...) 
whereas the risks of not treating it include the threat of blindness and major vessels occlusion. A patient 
with suspected GCA can receive a therapy based on Prednisone. Unfortunately, a treatment with high 
Prednisone doses may cause severe complications. Thus when confronted to a patient with suspected 
GCA, the clinician must adopt a clinical strategy. In |S] , the authors considered four different strategies: 

A : Treat none of the patients; 

B : Proceed to the biopsy and treat all the positive patients; 

C : Proceed to the biopsy and treat all the patients whatever their result; 

D : Treat all the patients. 

The clinician wants to adopt the strategy optimizing the patient outcomes measured in terms of utility. 
The reader is referred to m for more details on the concept of utility. The basic idea is that a patient 
with perfect health is assigned a utility of 1 and the expected utility of the other patients (not perfectly 
healthy) is calculated subtracting some “disutilities” from this perfect score of 1. These strategies are 
represented in Figures [4^ to [T^ with the different inputs involved in the computation of the utilities. 


A 


GCA 

9 


... 1 — dus — 

complication - . , 

ClUqc — aUdx 
gc 

No GCA 

complication - 1 — dug — dudx 

1 — gc 


No GCA 
l-i? 


1 — dudx 


Figure 4: The decision tree for the treat none alternative 


For example in strategy A (see Figure 4.31, the utility of a patient having GCA and developing severe 
GCA complications is given by 1 — dg — dugc — du^x- His entire sub-path is then 

g X gc X {I - ds - dugc - dudx)- 
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Prednisone 
p complication — 
Pc 


Test positive 
sens X g 


No Prednisone 
1— complication — 
1 -Pc 


GCA 

9pos 


No GCA 

1 - Spos 


GCA 

9pos 


No GCA 
1 - 9pos 


GCA 

gneg 

Test negative _ 

1 — sens X g 

No GCA 

1 9neg 


GCA 

complication 

gc{l - e) 

No GCA 
complication 
1 - gc(l - e) 


1 — dus — dup — 
dupc — dugc — du(, 

1 — dus — dup — 
dupc dtif^ 


1 — dup — dupc — dufy 


GCA 

— complication 

gc{l - e) 

No GCA 

— complication 
1 - jc{l - e) 


1 — dua — dup — 
dUgc — dUf, 


- 1 — dus — dup — duty 


1 — dup 


dui, 


GCA 

— complication — l-dus-dugc-dui, 

gc 

No GCA 

— complication — l - dug - du^ 

1-gc 


1 — dub 


Figure 5: The decision tree for the biopsy and the treat positive alternative 


Prednisone 
p complication 
Pc 


Test positive 
sens X g 


L No Prednisone 
complication 

1 -Pc 


Prednisone 
p complication 
Pc 


Test negative 
1 — sens X g 


No Prednisone 
1— complication 
1 - Pc 


GCA 

9pos 


No GCA 

1 - Spos 


GCA 

9pos 


No GCA 

1 - gpo. 


GCA 


No GCA 

1 gneg 


GCA 


No GCA 

1 gneg 


GCA 

p complication 
gc(l - e) 


1 — dug — dup — 
dXLpc dxigc dxib 


No GCA 
I— complication 
1 - gc{l - e) 


1 — dug — dup — 
dUpc — dub 


1 — dup — dupc — dub 


GCA 

— complication 

gc(l - e) 

No GCA 

— complication 

1 - Sc{l - e) 


1 — dug — dup — 
dUgc — dub 


- 1 — dug — dup — dub 


1 — dup 


dub 


GCA 

— complication 

pc(l - e) 

No GCA 

— complication 
1 - gc{l - e) 


1 — dug — dup — 
dUpc. — dUgr — dUh 


1 — dug — dup — 
dUpc — dUb 


1 — dup — dupc — dub 


GCA 

— complication 

gc(l - e) 

No GCA 

— complication 
1 - gc{l - e) 


1 — dug — dup — 
duy^ — dub 


- 1 — dug — dup — dub 


1 — dup 


dub 


Figure 6: The decision tree for the biopsy and the treat all alternative 
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GCA 

9 


GCA 

complication 
!)c(l - e) 


1 — dua — dup — 
dxtpc ” dUgc ” dUdx 


D 


Prednisone 
— complication 
Pc 


No Prednisone 
— complication — 
1 - Pc 


No GCA 
— complication 
1 - gc{l — e) 


1 — dUa — dUp — 
dup^ duif^x 


No GCA 
i-g 


1 - dudj: 


GCA 

9 


No GCA 
1-9 


GCA 

— complication 

gc{l - e) 

No GCA 

— complication 
1 -gc{l ~e) 


1 — dua — dup — 

dUgc - dudx 


- \ —dUa — dUp — dUdx 


1 — dup — dUfix 


Figure 7: The decision tree for the treat all alternative 


The input parameters 

As seen in Figures [T3| to [T3l the different strategies involve input parameters like e.g. the proportion g 
of patients having GCA or the probability gc for a patient to develop severe GCA complications (fixed 
at 0.8 as done in 0 ) or even the disutility associated to having GCA symptoms. Table summarizes the 
input parameters involved. 

The base values are provided by a physician expertise. The values P[-] and D{-) refer respectively to 
the probability of an event and to the disutility associated with an event. The minimum and maximum 
values m and M depict each parameter’s range for the sensitivity analysis. The base values are provided 
by a physician expertise. The utilities of the different strategies when all the input parameters are set to 
their base value are summarized in Table [21 

The base value of some input parameters are reliable while the others are really uncertain that leads us to 
consider them as random. As a consequence, if Ya, Yg, Yq and Yp represent the outcomes corresponding 
to the four different strategies A to D, the clinician aims to determine 

max{E [Y^], E[Yb] , E[Yc], E[Yo]} (18) 

with the uncertain model input presented in Table A sensitivity analysis is then performed to determine 
the most influent input variables on the outcome. 

Estimation phase and sensitivity analysis 

As done in [T2] and [5], all the random inputs will be independently Beta distributed. The Beta density 
parameters corresponding to each random input are determined by fitting the base value as their mean 
and capturing 95% of the probability mass in the range defined by the minimum and maximum. The 
remaining 5% will be equally distributed to either side of this range if possible. Concretely, each random 
input will be distributed as 

Z^m^Z<M + U)\rn>Z + VI Z^M 

where Z, U and V are independent random variables. Z is Beta distributed with parameters (a,/?). U 
and V are uniform random variables on [0, m] and [M, 1] respectively. 

Results 

The expected values of the utilities corresponding to the distributions given in Table[^are E[Y^] = 0.6991, 
E[%b] = 0.7570, E[Yc] = 0.7371 and E[Y£)] = 0.7171. Table summarizes the sensitivity measures of 
the seven random inputs with three different methodologies: considering the Sobol indices associated to 
the output vector Y = (Y^, Y^, Yc, Yp) (Multivariate) [T3] and the indices presented in this paper based 
on the Cramer von Mises distance. The last index considered is the one presented in |T| and named (3 
defined by 

Pi = E[sup{|Fy(j/) - FV|x.(j/)|}]. 
y(^y 
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Parameters 

Symbols 

Base 

Min. m 

Max. M 

Beta( 

a 

a,f3) 

P 

P[having GCA] 

9 

0.8 

- 

- 

- 


PJdeveloping severe complications of GCA] 

gc 

0.3 

0.05 

0.5 

4.179 

11.011 

P]developing severe iatrogenic side effects] 

pc 

0.2 

0.05 

0.5 

2.647 

10.589 

Efficacy of high dose Prednisone 

e 

0.9 

0.8 

1 

27.787 

3.087 

Sensitivity of temporal artery biopsy 

sens 

0.83 

0.6 

1 

7.554 

1.547 

D(major complication from GCA) 

dUgc 

0.8 

0.3 

0.9 

27.454 

6.864 

D(Prednisone therapy) 

dup 

0.08 

0.03 

0.2 

4.555 

52.380 

D(major iatrogenic side effect) 

dupQ 

0.3 

0.2 

0.9 

15.291 

35.680 

D(having symptoms of GCA) 

dus 

0.12 

- 

- 

- 

- 

D(having a temporal artery biopsy) 

dub 

0.005 

- 

- 

- 

- 

D(not knowing the true diagnosis) 

dudx 

0.025 

- 

- 

- 

- 


Table 1: The data used by Buchbinder and Detsky [8] in their analysis 


Treatment alternative 

Utilility 

A Treat none 

0.6870 

B Biopsy and treat positive 

0.7575 

C Biopsy and treat all 

0.7398 

D Treat all 

0.7198 


Table 2: 


The utilities of the different strategies when all the input parameters are set to their base value 


Sensitivity meas. 

Ranking 

Multivariate 

1236475 

Baucells et al. 

1627354 

Cramer von Mises 

1627354 


Table 3: Sensitivity measures 
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We then use the estimator given in O Table 1] adapted to the multivariate case that is based on the 
tedious and costly estimation of conditional expectations. 

As a conclusion, both methodologies based on the whole distribution provide the same ranking unlike the 
multivariate sensitivity indices. Nevertheless, the main advantage of the Cramer von Mises sensitivity 
methodology is that one can use the Pick and Freeze estimation scheme which provides an accurate 
estimation simple to implement. In |S], the authors study a slightly different model that explains the 
numerical differences between their results and the ones of the present paper. Furthermore, they perform 
a sensitivity analysis on the best alternative with the greater mean instead of considering the multivariate 
output. 
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